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Abstract 

The dynamics of the dissipative two-level system at zero temperature is studied using three 
different cumulant expansion techniques. The relative merits and drawbacks of each technique 
are discussed. It is found that a new technique, the non-crossing cumulant expansion, appears 
to embody the virtues of the more standard cumulant methods. 
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1 Introduction 



The standard spin-boson problem, described by the Hamiltonian 



H 



A 
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has served as a paradigm for the description of dissipative effects in condensed phases. Some 
experimental realizations of such a Hamiltonian include, e.g., the detection of macroscopic quantum 
coherence in superconducting quantum interference devices |3j, |3| , tunneling effects in metallic and 
insulating glasses M, electron transfer reactions || and the diffusion of light interstitial particles 
in metals J?]]. In each situation, the physical realization of the parameters in the Hamiltonian (1) 
is different. For instance, in metallic glasses at low temperatures the electron- hole pairs at the 
Fermi level constitute the bosonic bath, while for insulating glasses, tunneling effects are damped 
by localized and delocalized vibrational modes. Thus, the Hamiltonian (1) embodies a wealth of 
physical situations and has been studied in great detail (see for instance [|], |2| and references quoted 
therein, or more recently ||, ^0|). 

In order to study the dynamics of the two-level system coupled to a harmonic bath as in (1), we 
need a method of "tracing out" the bath or spin degrees of freedom. The bath degrees of freedom 
can be specified by the spectral density function, 



rrijUj 



-5{ljJ — LOj), 



(2) 



which gives the bath density of states weighted by the square of the coupling strength between the 
two-level system and the bath. In most studies of the spin-boson problem, the spectral density 
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takes the ohmic form [|], §], 

J(u) = 2irau exp(— u/u c ), (3) 

where a is a measure of the coupling strength, and ui c is a frequency cutoff for the bath. We note 
that in many cases, such as the coupling of a spin degree of freedom to a three dimensional phonon 
bath in the deformation potential approximation, the spectral density (3) is not realistic, and must 
involve higher powers of u. 

The usual approach to finding reduced equations for the spin variables of interest involves the 
use of the functional integral formulation of quantum dynamics [|l], Q. Formally exact equations 
may be found for the variables 

P(t) = Mt)), (4) 

and 

C(t) = \{M0),a M {t)}) Pl (5) 

where (•••) ( g refers to an average with respect to the canonical ensemble of (p]). The quantity P (t) 
describes the population difference in the localized spin states of the Hamiltonian (1), given that 
the particle is initially localized in one well and in thermal equilibrium with the bath. It is the 
variable of interest in certain physical situations, for example, the electron transfer problem ||. 
The quantity C(t), the symmetrized equilibrium correlation function of the tunneling coordinate, 
is related to the structure factor for neutron scattering off the tunneling particle, and is of great 



significance in various problems, including the antiferromagnetic Kondo problem [11|. For C(t) the 
long-time behavior at zero temperature is known from the generalized Shiba relation which predicts 
algebraic decay C(t) oc t~ 2 |2], [l^]. For P(t) the situation is less clear, however some studies have 



predicted exponential decay as t —* oo 0, 13 1. Despite the importance of C(t), we will focus on the 
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variable P{t) in the following. 

The formal path integral expression for P(t) is extremely cumbersome, and a suitable approx- 
imation must be implemented to obtain useful information. The so-called "non-interacting blip 
approximation", or NIBA (2|, is the most commonly used approximation. In this scheme C(t) is 
entirely determined by P(t), i.e., C(t) = P(t). The NIBA may be obtained from the exact expres- 
sion for P(t) by invoking a series of physically based approximations. For very low temperatures, 
these approximations often break down, unless a is very small and only short times are considered. 
At zero temperature, the NIBA is not justified in the antiferromagnetic Kondo regime \ < a < 1. 
The NIBA also incorrectly predicts asymptotically algebraic, rather than exponential, decay for 
the variable P(t). Lastly, NIBA incorrectly predicts that at zero temperature, C(t) ~ t^ 2 ^^ a \ 

Despite these flaws, the NIBA is useful for obtaining quantitative results for P(t) for high 
temperatures, when the tunneling dynamics is incoherent, and in predicting the qualitative behavior 
of P(t) for low temperatures. For instance, at zero temperature, the NIBA correctly predicts a 
crossover from damped oscillations to incoherent decay for the variable P(t) at the point a = i. 



As shown by Aslangul et. al. [14], the NIBA may be obtained by first applying a small po- 
laron transformation to the Hamiltonian (1), followed by a second order application of the usual 
Nakijima-Zwanzig projection operator technique. It has been known for some time that this pro- 
jection technique, which leads to a master equation of the convolution form, is an order by order 
resummation of a particular type of cumulants known as "chronologically ordered" cumulants 
]i~6| , [l7], [l8|, 19, 2C]. The use of the "chronological ordering prescription", or COP, when truncated 
at second order thus leads to the NIBA. 

Interestingly, Aslangul et. al. [^] earlier applied a convolutionless master equation technique 
to the study of the zero temperature spin-boson problem. This type of master equation, which can 
be derived by using a different type of projection operator, involves the summation of a different 



type of cumulants, known as "partially (time)-ordered" cumulants [21, 17, 18|, 19, 20]. This method 
was probably abandoned for two reasons. First, it incorrectly describes incoherent relaxation for 
P{t) for all values of a. Secondly, it cannot be obtained in a simple manner from the exact 
path integral expression. The second objection is irrelevant, since it is still possible that such an 
approximate resummation describes the exact behavior of P{t) well. The first flaw, however, is 
quite serious. Despite this, the expression obtained from the "partial ordering prescription", or 
POP, which naturally resums to an exponential form, may be expected to give a better description 
of P(t) in the incoherent region. In fact, for values of a greater than |, but not too large, this 
method indeed describes (nearly) exponential relaxation. Furthermore, as will be demonstrated in 



this paper, recent simulations of Egger and Mak [22| show that the POP method more accurately 
captures the deep decay of P(t) at zero temperature for a > ^ than does the COP (NIBA) method, 
even before the algebraic behavior of the NIBA is manifested. 

It is well known that by choosing a particular ordering in a truncated cumulant expansion, we are 
implicitly assuming different statistical properties for the relevant bath operators. The first purpose 
of this paper is to specify these statistical properties for the case of the spin-boson problem at zero 
temperature. Using this "stochastic" type intuition, we then discuss various cumulant ordering 
schemes and their associated descriptions of the behavior of P(t) at T = 0. This paper is organized 
as follows: In Sec. 2 we first present a new derivation of the exact expression for P(t) that allows 
for clear specification of the statistical properties of the bath. For this purpose, orthogonally to 
the conventional approach, we first integrate out the spin degrees of freedom exactly. In Sec. 3, we 
briefly discuss the COP and POP methods. We then turn to a recently introduced new cumulant 



method, the "non-crossing" cumulants [23, 24, 25|. Lastly, in Sec. 4, we compare the methods to 
exact simulation result. 
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2 Moment Expansion 

We begin with an explicit expression for P(t) through fourth order in A. We could, if we wished, 
obtain these terms from the exact path integral expression for P(t), however, we believe that the 
method used in this section most clearly shows the connection to the stochastic methods upon 
which the cumulant expansions are based. In effect, our method offers another route to the formal 
expression of Ref. ||, |2| . 

We begin with the Hamiltonian (1) in the form, 



H = H' + ^Px, (6) 

9l 



k ^ 



The quantity we wish to calculate is P(t) which is defined as 



P(t) = (a z (t)) = Z~ lr Fr (exp(iHt)a z (0) exp(-im)vr + exp(-/3#V + ) (8) 



where 



Z = Tr (7r + exp(-/3i7')vr H 
a z = \L)(L\-\R)(R\, 
a x = \L)(R\ + \R)(L\, 



and (5 is the inverse temperature. We now diagonalize (6) in the spin manifold with the use of 
a transformation employed by Shore and Sander [26, 27] in their study of the self-trapping of an 
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exciton coupled to phonons, namely, 



U 



V2 



(9) 



where 



(-l)E k b W =exp (m^blh 
\ k j 



The operator (f> is seen to be the parity operator for the bath modes. In the transformed picture, 
we can express 



P(t) = -Z _1 Tr (exp(iHt)a x {0) exp(-iHt)n + exp{-l3H')n + ^j 



(10) 



where 



H 



— / 

H 



A 



>cr z + H , 



71 = o( 1 - C7 ^)> 



(11) 

(12) 
(13) 



and Z is now defined with respect to H and ir + . 

We now perform the trace over the spin degrees of freedom in (10), leaving 



P(t)=Re[G(t)}, 



(14) 



where 



G(t) = Tr b (exp(iH+t) exp(-iif-t) exp(-/?tf')) / Tr 5 (exp(-/?F')), (15) 
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with 

A 

H ± = ±-<j> + H. 

This trace over the bath degrees of freedom is most easily performed in the small polaron repre- 
sentation, defined by the transformation 



U = exp(0, (16) 

k ^ 



In this picture, we may express G(t) as 



(exp^(i / dTi]{r)) exp^(i / (Iti](t))}b, (18) 

JO JO 



where 



r ? (t) = ^expK(t)]<Aexp[^)], (19) 

and 

£(t) = V ^L(6 fce - <w ** - fete*"**). (20) 

The averaging (denoted by (...)s) is over the canonical ensemble of harmonic oscillators, (i.e., pb = 
exp(— (3 J2k Wkb^bk) / Trj J exp(—(3'^2 k uJkb k bk), and exp^ (exp^) denotes a time ordered exponential 
with latest time to the right (left). From this point on, all averaging will be with respect to this 
ensemble, and we will drop the subscript "B" . Since the spin degree of freedom has been removed, 
our method allows us to focus on the bath operators that arise in the expansion of P(t). Using the 
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following properties of the parity operator, 



exp[£(t)] = exp[-£(t)] . 



and 



6 2 = 1, 



we can show, through fourth order in A, the moment expansion for P(t), 



pt /"^i /*^2 /*^3 

P(t) = l+ dh dt 2 m 2 (t 1 t 2 ) + dh dt 2 dt 3 dt 4 m 4 (h,t 2 ,t 3 ,t 4 ) + (21) 
Jo Jo Jo Jo Jo Jo 



where the moments rrii equal 



m 2 (h,t 2 ) = -A 2 Re(B_(ti)S+(t 2 )), (22) 
A 4 r 

m4(h,t 2 ,t 3 ,U) = — Re [(B_(ti)5+(t 2 )B_(t 3 )5 + (t4)) 
+ (B_(t 2 )5 + (ti) J B_(t 3 )fl + (t 4 )) 
+ <S_(t 3 )B + (ti)S-(*2)B + (t 4 )> 

+ (B_ (t4)B + (ti)B_(t 2 )B + (t 3 )) ] , (23) 



and 

B±(t) = exp[±2£(*)]. (24) 

Note that m-2 n -i = 0. In this paper, we shall only use the first two nonvanishing moments, although 
it is a simple matter to execute the expansion to an arbitrary order. From (^l|)-(|23|) we conclude 
that P{t) is entirely determined by the statistical properties of the bath operators B±(t) with 
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respect to the canonical state of the bath. Note that the operators B± always appear in pairs. In 
order to specify the statistics obeyed by the operators B±, we now calculate the second and fourth 
moment of the B± , s. It is a simple matter to show that 



(B_(*i)B+(* 2 )) = exp HQi(ti - t 2 ) - Q 2 (h - t 2 )} , (25) 



where 



Qi(h - t 2 ) = 4£ ( ^-Y erin(w fc (ti - t 2 )), (26) 



A' 



and 

Q 2 ( h - t 2 ) = 4^2 ( —) [1 - cos(w fc (ti - t 2 ))] cothC9w fc /2). (27) 

Furthermore, by using the relation e A+B = e A e B e~^ A > B ^ and (e Kbt e~ K * 6 ) = e~' K l 2 for [A, JB], k being 
c-numbers, one can also show that the operators B± have the following "statistical" property, 

(fl-Ct!)^^.^^)) = (B^B^MB.it^iU)) • 

(28) 

This property can be extended to an arbitrary number of B± pairs. This gives a type of "Wick" 
theorem for the operators B±, and demonstrates the underlying reason why only the functions Q\ 
and Q 2 appear in the exact path integral (see Eqs. (4-17) to (4-22) in Ref. M). It can now be 



explicitly checked that the expression (21) is identical to the exact path integral expression, at least 



through the fourth moment. Note that the property (28) is different than the statistical properties 
held by commonly used stochastic processes such as Gaussian, two-state-jump, or Gaussian random 
matrix processes. We will return to this point in the next section. 

The moment expansion itself is not a very useful scheme for describing dynamics, because an 
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arbitrary truncation of the expansion leads to secular terms that grow with time. We next resort to 
schemes that provide partial (approximate) resummations of the moment expansion to all orders. 
Such schemes are the cumulant expansions that will be introduced in the next section. 



3 Cumulant Expansions 

We now discuss the various ordering prescriptions which allow for partial resummation of the 



expansion (21). Each ordering method leads to a unique type of master equation [24]. We note 
that, when carried out to infinite order, all of the ordering techniques give the same (exact) result. 
When truncated at a finite order, however, the results are different. In simple stochastic situations, 
when the generator for time evolution (the Liouville operator) commutes with itself for all times, 
i.e., [L(t),L(t )] = 0, the use of a particular truncated cumulant expansion implies a knowledge 
of the stochastic properties of the bath functions. In simple cases, truncation of the cumulant 
expansion in the "correct" ordering prescription can lead to exact results that may be obtained in 
the "incorrect" ordering prescription only at infinite order. In the quantum case described by the 
Hamiltonian (1), where [L(t),L(t')] ^ 0, truncation of a cumulant expansion at finite order in any 
ordering prescription will never lead to exact results due to the non-commutivity of the Liouvillian 
at different times [^, ^Cj] . It is precisely this noncommutivity that leads to the variety of 
time-orderings of the operators B± in the expression (p3|) for 7714. Despite this fact, the statistical 
properties of the bath operators still dictate the choice of an ordering prescription that provides 
the most rapid convergence of the cumulant series (if such convergence exists) plj . 



We begin by discussing the chronological ordering prescription, or COP. In this prescription, a 
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master equation of the form (see for instance [2C 

dP(t) 



dt 



[ K COP (t,T)P( T )dT (29) 
Jo 



is obtained. The kernel K COP (t,r) is obtained from the moment expansion as 



K COP (t,T)=j 2 (t,T) + Y / dr 2 ... dr nln+1 (t,n,...T n ). (30) 

n=2 Jo Jo 



The COP cumulants, 7, are obtained from the moments by a recusion relation [24|. In the present 
case this yields 

72„-i = 0, 
l2(t,n) = m 2 (t,ri), 
74(i,Ti,T2,r 3 ) = m 4 (t,Ti,T2,T 3 ) - m 2 (i, n) m 2 (r 2 , t 3 ) , . (31) 

For the simple case where the stochastic Liouvillian commutes with itself for all times, all of 
the COP cumulants 7„ vanish for n > 3 if the stochastic bath functions have the two-state-jump 
behavior [^0] 

{B(t{)B(t 2 )B(h)B{t A ) ■ ••) = {B{t l )B(t 2 )){Bm)B{U) ■ ••), (32) 

for t\ > i 2 > £3 > £4-", where -B(i) is the stochastic bath function responsible for system dissipation. 
If these bath functions have different statistics, it may not be a good approximation to truncate 
the series at low orders. 

Returning now to the quantum case of interest in this paper, we find at lowest order, as shown 
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by Aslangul et. al. @, the NIBA equation for P(t) 

^ = J*m 2 (t-T)P(r)dT, (33) 
where, at T = 0, using the ohmic constraint (|3|), along with (p2|), we may express 

m2 ( t - T) = -A^ T — 1_ . (34) 



As shown by Grabert and Weiss [p^ l, the solution to ( p3| ) with the kernel ([34|) can be given for all 
a < 1 (in the limit ^ — > 0) by the Mittag-Leffler function |^| , 

P™(y) = ^ 2( i- a) (-y 2(1 - a) ), (35) 

where y = A e ffi and 

A cff = A [cos(vra)r(l - 2a)] ^a^T (— \ 1_a . (36) 



This solution shows damped oscillations for a < i, and incoherent decay for a > i. This behavior 



has been qualitatively confirmed by Monte Carlo simulation [22]. As mentioned in the introduction, 
the NIBA cannot give the correct asymptotic decay of P(t) (yielding the algebraic decay P(t) oc 
£-2(i-a) pathgj- than exponential decay), and is unable to account for the depth of the decay in the 
region a > \ even before the incorrect algebraic behavior sets in. The NIBA is, however, known 
to work quite well for short times and weak coupling strengths. The analysis given in the last 
section provides a novel explanation for this fact. For "short" times and "small" values of a the 
function rri2{t) is a rather broad, weakly decaying function of time. When this is the case, the 



statistical property (|28| ) of the operators B± is approximately of the two-state-jump form (32) as 



12 



far as the integrations over the cumulants r y n >3 ar c concerned. This approximate equivalence holds 
in a stochastic sense, in that all of the four point correlation functions in m\ (see Eq. (p3)) may 
be approximated by 77t2(£i, t-i) m-iit-A, £4)- For such times and coupling strengths, the NIBA will be 
essentially exact, as all COP cumulants for n > 3 will vanish when integrated. We shall not provide 
precise meaning to the terms "short" or "small", although their meaning should be clear in the 
context of the present discussion, and could be quantified without undue labor (in fact "short" and 
"small" will be coupled in the sense that the effective timescale of oscillation or decay, (A c ff) _1 , 
depends on a). Note that the statistical property ( p8[ ) trivially gives two-state-jump behavior for 
a = 0, which leads to the correct behavior P(t) = cos(At). While this is obvious, other cumulant 
techniques (such as those discussed below) do not embody this type of statistics for a = 0, and 
cannot give the correct, freely oscillating solution for zero coupling strength upon truncation at 
second order. The statement (often given in the literature |34j| ) that NIBA works for weak coupling 
because it is a perturbative scheme is thus not strictly correct. 

The (somewhat heuristically) demonstrated fact that the property ( |28| ) can resemble two-state- 
jump behavior under certain circumstances leads one to believe that extending the COP scheme 
to fourth order would not be useful, since this property is reflected in the vanishing of all COP 
cumulants higher than the second. Extending the COP method to fourth order does not give a 
method for computing "interblip " interactions in the language of Ref . |l]] . 

We now turn to the partial ordering prescription, or POP. At second order, this method was 
applied by Aslangul et. al. [15] to the spin-boson problem at T = 0. The POP master equation 



has a convolutionless 24] form 



dP(t) 
dt 



J* K POP { T )dT^j P{t). (37) 
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K POP (t) may be obtained from the moments 

K POP (t) = J2 M T i / dr 2- / dT n 6 n+1 (t,n,...T n ), (38) 

~L JO JO JO 

where 

#2n-l = 0, 
02(t,n) = m 2 (t,Ti), 
^4(*,ti 3 7"2 3 t 3 ) = m 4 (t,ri,T2,r 3 ) - m 2 (t,ri)m2(r2,T 3 ) 

- m2(t,r 2 )m 2 (ri,T3) - m 2 (t, r 3 ) m 2 (ri, r 2 ), (39) 

The POP resummation is exact at second order for the simple case of a classical Gaussian 
stochastic process. We note that the statistical property (^q) appears to be very different from the 
standard Wick theorem for Gaussian processes. We may still expect that the POP method is better 
suited for the incoherent regime a > ^ for the following reasons. First, the POP technique resums 
to an exponential form, which is expected to better capture the long time behavior of P(t). In 
general, the POP method sums (infinitely) more terms than the COP method does. For example, 
expansion of the second order truncation in the COP gives, to fourth order 

ft rt± rt rti rt 2 rt-i 

P(t) = l+ dh dt 2 m 2 (h,h) + / dh / dt 2 dt 3 dUm 2 (t 1 ,t 2 ) m 2 {h,U) + - 
J Jo Jo Jo Jo Jo 

whereas the POP gives 

P(t) = 1+/ dh dt 2 m 2 [t 1 ,t 2 )+ I dh dt 2 dt 3 dt 4 m 2 {h, t 2 ) m 2 (t 3 , U) 
Jo Jo Jo Jo Jo Jo 

+ rn 2 (h,t 3 ) m 2 (t 2 ,U) + m 2 (h,U) m 2 (t 2 ,t 3 ) + ... . 

14 



Clearly, the extra terms do not insure a more accurate result. For example, for weak coupling 
strengths, the POP must be carried out to infinite order to obtain coherent behavior. However, in 
the incoherent regime, the effective timescale, defined by (|36| ) is very long, while the decay of the 
function 777.2 (i) is "slow" (algebraic). In this case, we may expect that we must include terms like 
^2(^15^4)^2(^2^3) that extend over large portions of the integration region. As we will show in 
the next section, the POP method captures the behavior of P(t) better than the COP method in 
the incoherent regime, even before the full asymptotic time behavior is displayed. 

At second order, (K POF '(t,r) = miit — t)) the POP equation (38) may be solved [15], 



Pit) = exp 



A 2 1 / cos(2(l - a) tan' 1 Lo c tY 

4u>* (a-l/2)(l-a) 1 1 (l+c^ 2 ) " 1 



(40) 



Note that equation fl40| ) describes a stretched exponential rather than exponential decay. For values 
of a that are not too much larger than ^, however, the POP expression should be much more 
accurate than the COP expression, at least at second order. 

We have now given some motivation for the belief that the COP method (at lowest order) should 
give a better description of P(t) in the region a < \ while the POP method should be better in 
the incoherent region a > | • We now ask whether there is a summation method that is a "hybrid" 
of the two methods, in the sense that it can incorporate at low order features of the COP and 
POP methods. In the theory of stochastic processes, such a technique has recently been developed 



[23, 24, |2q|. This method is based on the summation of "non-crossing" (NC) cumulants (for a 



precise definition see Refs. [23, 24, ^5|). For simple stochastic situations, if the coupling is not 
too strong, the NC technique (including terms up to fourth order) has been shown to interpolate 
between the two-state-jump behavior and the Gaussian behavior 24]. 



The NC description leads to a nonlinear equation of motion for Pit) ]24|1, which at second order, 



15 



may be expressed 

dP(t) 
dt 

where 



M(t), (41) 



M(t) = dh C 2 (t - ti)P(t - *i)P(ti). (42) 

JO 

To fourth order the master equation for P(t) in the NC scheme reads 
dP(t) rt rti rta 



dt 



M(t) + { dh f 1 dt 2 [' dt 3 ( 4 (t, t 1 ,t 2 ,t 3 )P{t - tx)P(h - t 2 )P(t 2 - t 3 )P(t 3 ). (43) 
jo Jo Jo 



As in the previous two case, the NC cumulants, £, may be obtained from the moments by a recursion 



relation [23, g^]. In the present case this yield up to fourth order 



C2n-1 = 0, 

C 2 (Mi) = rn 2 (t,ti), 
(4{t,h,t 2 ,t 3 ) = m4(t,t 1 ,t 2 ,t 3 )-m 2 (t,t 1 )m 2 (t 2 ,t 3 )-m 2 (t,t 3 )m 2 {t 1 ,t 2 ),... . (44) 

It is clear that in appearance, the NC cumulants are a "compromise" between the COP and 
POP cumulants. We note in passing two interesting facts. First, in the stochastic realm, the NC 
ordering prescription truncated at second order is exact for the case of a stochastic bath modeled by 
symmetric (N x N) Gaussian random matrices for the P's. In this case, the "crossing contraction" 
m2(t, t 2 ) m 2 (ti, t 3 ) vanishes by means of a i argument for N — > oo. This leads naturally to the 
equation (^) and its systematic generalization (|43| ) through the NC cumulants. Eq. (f42|) has first 



been derived by Kraichnan 36, 37, g8|, 3£] in the field of turbulences and fluid dynamics. Our 



motivation for the application of this method is not based on a stochastic type of reasoning, but on 
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the fact that in simple situations this ordering prescription may combine the benefits of the COP 
and POP methods. 

Before concluding this section, we would like to apply all three ordering prescriptions to the 
case T = 0, a = ^. Here, it is known that the "exact" (in the sense specified in Ref. Q) result for 
P(t) is 

p (( )= CX p(-M£0 (45) 

in the limit £j- — ► 0. Note that in this limit, the second moment becomes 5-correlated (t > 0) 

Using the fact that P(0) = 1, it is clear that all three ordering prescriptions give the same result 
given by fl45|) at second order. Hence, the value a = \ corresponds to the white noise limit of the 
bath operators B±(t). 

4 Results and Conclusions 

Before comparing the results of the three ordering methods, we make some comments on the 
methods discussed in Sec. 3. We have shown how three different cumulant methods give rise to 
different master equations with different properties. We have tried to physically motivate when 
each approach should have success when applied to the spin-boson problem at zero temperature. 
Note that in general, the discussion of convergence of each cumulant series is a difficult task. This 
task is made more difficult by the fact that, at zero temperature, the algebraic decays of the bath 
correlation functions leave us with no clearly defined relaxation time for the bath. This means that 
we will rely almost exclusively on physical considerations and comparison with accepted results 
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to determine the success or failure of the methods employed. The case of finite temperature, 
which can be studied by the same methods employed here, is often easier in this respect. If an 
exponential correlation time t\, can be assigned to the decaying bath correlation functions, than it 
is possible to consider a systematic expansion in A e fjTb provided that this dimensionless parameter 
is small. When this is the case, the POP provides the most facile way of systematically summing 
terms in the parameter A^tj, |HJ]. In case of ohmic dissipation and finite temperature T, the 
characteristic correlation time of the B±(t) is given by t& ~ (2-KaT)^ 1 This point of view 
provides a novel explanation for the familiar statement that the NIBA works well in the incoherent 
tunneling regime A e ffTfr < 1, In a stochastic language, this parameter region corresponds to the 
narrowing or Markov limit of the -B±'s. Similarly to the white noise limit mentioned above, one 
finds that all three cumulants schemes work well already at second order and provide essentially 
the same behavior, P(t) « exp[-A 2 eS (T)n] with A cfT (T) oc T a Q. 

Since we expect the NIBA to be accurate for very weak coupling strengths, we first turn to 
the case of weak to intermediate coupling strength, a = 0.3. For coupling strengths in this range, 
simulations at low temperatures have shown that the NIBA is qualitatively correct in predicting 
damped oscillations, but may fail in predicting the damping strength. An example of this is given 



by the simulations of Makarov and Makri [13] which show that for intermediate coupling, NIBA 
may fail by slightly underestimating the number of oscillations in P(t). We note, however, that 
these simulations were carried out for values of that are not very small. In Fig.l, we plot the 
NIBA (second order COP) solution for P(t) against the solutions obtained from second and fourth 
order truncations of the non-crossing cumulant method, and the second order POP. Note that, 
as expected, the second order POP solution for P(t) fails to produce any oscillations. We expect 
that for a < \ the POP will always be inaccurate at low orders. The second order non-crossing 



cumulant solution for P(t), obtained from the Kraichnan-type equation fl4l|)-([42|) is similar to the 
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NIBA solution, although the oscillation in P(t) is much weaker. The fourth order noncrossing 
cumulant solution gives a first oscillation which is very similar in magnitude to the NIBA solution, 
however, it describes one extra weak remnant of an oscillation. This behavior is very similar to 



the behavior displayed in the exact simulations of Makarov and Makri 13 ] . Although this example 
represents only one value of a, similar results may be obtained for all moderately strong values 
of a up to a = ^. Thus, it appears that the non-crossing scheme works well in incorporating the 
qualities of the COP method for moderate values of a when a < 5. 

One may be a bit alarmed by the magnitude of the difference between the second and fourth 
order plots of P(t) obtained via the non-crossing cumulant method. This is not necessarily indicative 
of a lack of convergence in the summation of the cumulants. As an example to support this, we 
mention the cumulant expansion results of Aihara, Budimir, and Skinner |3l| . For a different 
model of relaxation, these authors compared the results obtained from truncation of the POP 
expansion at sixth order with exact simulations. For coupling strength that were not too large, 
they found that by sixth order, the cumulant method quantitatively describes the results obtained 
from simulations. In cases where convergence was obtained, there was a large difference between 
the solutions obtained from truncating the cumulant series at second and fourth orders. In such 
cases, the fourth order solution was nearly quantitative, but sightly overestimated the exact result. 
We believe that this is precisely what is occurring here, although we have no direct evidence (i.e., a 
sixth order result) for this belief. We base this statement on the similarity with exact simulations 
at slightly different coupling strengths, the fact that the fourth order non-crossing result shifts P (t) 
in the right direction, and the reasonable appearance of the result. 

We now turn to the relaxation of P(t) in the incoherent regime a > i. Here, the beautiful path 
integral simulations of Egger and Mak [^] provide a means of comparing the cumulant expansion 
methods with exact results. In this region, we expect the POP to be most successful, while the 
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NIBA (second order COP) is expected to be worse. Based on experience with simple stochastic 
situations, we hope, as in the coherent portion of the coupling space, that the non-crossing scheme 
can capture the essence of the POP in this regime. In Fig. 2, we show the decay of P(t) calculated 
by differing ordering prescriptions for a = 0.6. It must be noted that the simulations were carried 
out for long times, but not long enough to show the asymptotic algebraic decay of the NIBA (second 
order COP) solution of P(t), or the asymptotic exponential decay of the exact solution. Regardless, 
the second order truncation of the POP still gives the most accurate description of the decay of 
P(t). As we had hoped, the second order non-crossing technique is nearly identical to the POP for 
this value of a. Figure 3 shows the results for a = 0.7. Again, the POP performs the best, while the 
second order non-crossing scheme over estimates the decay. As in the case of weaker coupling, we 
see if truncation after fourth order in the non-crossing cumulants can properly correct the second 
order result. This test is shown in Fig. 4. While the results appear to show that the non-crossing 
scheme is converging to a POP-like description of P (i), we again must exercise caution due to the 
lack of further concrete evidence for this belief. For such large value of a, it is quite possible that 
the cumulant methods break down. 

One interesting property displayed in Fig. 3 and Fig. 4 is the near quantitative agree- 
ment between the POP description of P(t) and the exact simulation of P{t) for moderately long 
times. In order to test if this is a coincidence, we have computed dy Hf dy K POP {y ) and 
Jq dy Jq dy K POP (y ) where K POP (t) is the n-th term in the expansion fl38|). If the integrated 
second order POP cumulant is of order one for a given time interval, while all other POP cumu- 
lants are small when integrated over the appropriate time domain, then we expect the truncation 
at second order to be a good approximation. While we cannot study all the POP cumulants, we 
have studied the second and the fourth. In Fig. 5, we compare the properties of the second and 
fourth POP cumulants for a = 0.7. For y = 1.4 to y = 2 (the boundary of the simulation results 
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of Egger and Mak |22|| ), we see that the contribution from the second POP cumulant is at least 
ten times greater than the contribution from the fourth cumulant. This strongly suggests that the 
agreement of the second order POP method with the exact simulations is no coincidence. In fact, 
the agreement between (41) and the simulation occurs precisely in the interval where the second or- 
der cumulant dominates the fourth order cumulant. Since the slopes of the two curves suggest that 
this behavior continues for some time, we feel there is strong evidence for the somewhat remarkable 
conclusion that, for significant intermediate times, the decay of P(t) is quantitatively described by 
a stretched exponential. For longer times, the decay is most likely purely exponential. 

We now summarize the results presented in this paper. We first carried out a new derivation of 
the moment expansion for the variable Pit) in the spin-boson problem. We then used this derivation 
to discuss the "statistical" properties of the relevant bath operators. Using the moment expansion, 
we first discussed the chronological and partial ordering prescriptions that involve different types 
of cumulants. We discussed the merits and drawbacks of each method. In an effort to combine 
the merits of the COP and the POP, we applied the non-crossing scheme. Specializing to the case 
of zero temperature, we tested each method, including fourth order terms when necessary. Our 
results show that the non-crossing scheme is a promising candidate for combining the virtues of 
the COP and POP, especially for intermediate values of a on either side of the coherent-incoherent 
transition value of a = ^. We note that more work should be done to test the validity of this claim. 
Lastly, we have provided evidence to support the belief that the stretched exponential behavior 
described by second order truncation of the POP in the incoherent portion of coupling space may 
infact be very accurate for intermediate times. 
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Figure Captions 

Fig.l: Zero temperature plot of P (y) (y = A c gt) for a = 0.3 and ^ = 6. The dotted line is the 
second order POP result, the dashed line is the second order non-crossing cumulant result, the 
dash-dotted line is the NIBA (second order COP) result, and the solid line is the fourth order 
non-crossing cumulant result. 

Fig.2: Zero temperature plot of P(y) (y = A c fjt) for a = 0.6 and ^ = 6. The dash-dotted line is 
the NIBA (second order COP) result, the dashed line is the second order POP result, the solid line 
is the second order non-crossing cumulant result, and the open circles are the simulation result of 
Egger and Mak f22j. 

Fig.3: Zero temperature plot of P(y) (y = A c fjt) for a = 0.7 and ^ = 6. The dash-dotted line is 
the NIBA (second order COP) result, the dotted line is the second order POP result, the solid line 
is the second order non-crossing cumulant result, and the open circles are the simulation result of 
Egger and Mak f22j. 

Fig.4: Zero temperature plot of P{y) (y = A c fji) for a = 0.7 and ^ = 6. The dotted line is the 
fourth order non-crossing cumulant result, the dashed line is the second order POP result, and the 
solid line is the second order non-crossing cumulant result. Note the change in the x-axis. 
Fig. 5: Relative magnitude of second and fourth cumulant effects in the POP for a = 0.7. The 
dashed line shows \ Jq dy Jq dy K\ (y )| and the solid line shows \ Jq dy Jq dy K^ OP (y )|. 
K POP (t) is defined in Eq. (38). 
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